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Abstract: We present a new multi-channel integration method and its implemen- 
tation in the multi-purpose event generator MadEvent, which is based on MadGraph. 
Given a process, MadGraph automatically identifies all the relevant subprocesses, 
generates both the amplitudes and the mappings needed for an efficient integration 
over the phase space, and passes them to MadEvent. As a result, a process-specific, 
stand-alone code is produced that allows the user to calculate cross sections and pro- 
duce unweighted events in a standard output format. Several examples are given for 
processes that are relevant for physics studies at present and forthcoming colliders. 
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1. Introduction 

The high energy and luminosity of the present and future colliders, from the Tevatron 
to the Large Hadron Collider to the proposed 500 GeV Linear Collider, will offer the 
widest range of physics opportunities for the exploration of the high-energy frontier. 
Among the highest priorities is understanding not only the nature of the electroweak 
symmetry breaking but also the mechanism through which the electroweak scale sta- 
bilizes. Simple and very robust arguments indicate that this scale should be less than 
one TeV, very much in the reach of the above mentioned colliders. At this energy, 
multijet final states in association with leptons or missing energy constitutes the most 
interesting data samples. The quest for the Higgs boson(s) and/or for supersymmet- 
ric particles will rely on our ability to predict and understand the standard-model 
processes which are the backgrounds to these searches. In general, not only pure 
QCD interactions lead to multijet configurations, but also the production of vector 
bosons or top quarks that subsequently decay into quarks yield similar signatures 
with many particles in the final state. 

Although straightforward in principle, the evaluation of the cross sections and 
especially the generation of unweighted events for such processes is not a trivial task. 
First of all, the number of subprocesses contributing to a process can be fairly large. 
For instance, for W + 4-jet production in pp collisions there are 498 distinct parton- 
parton subprocesses. Clearly, it is desirable to automate not only the identification 
of these subprocesses, but also their integration and event generation. 

Unfortunately, the complexity of each subprocess also grows quickly with the 
number of particles in the final state. If the amplitudes involve identical particles then 
the growth is factorial. A typical example, amply discussed in the literature [1], |], 3]) 
is the computation of n-gluon amplitudes where already for a relatively small number 
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of gluons produced, say n > 6, clever reorganizations of the perturbative expansion 
have to be employed to reduce the factorial growth to a more manageable exponential 
one[|,§. 

The second step is the integration of the squared amplitude over the phase space, 
which is necessarily performed using Monte Carlo techniques. This task involves 
solving several issues. In the first place, the amplitude results in a very complicated 
function of the momenta which displays sharp peaks in many region of the phase 
space. In order to achieve an efficient numerical integration, one needs to identify 
the position of these peaks and to map them onto several different sets of variables, 
which are usually referred to as "channels" . In principle, both tasks are not trivial. 
Fortunately, a standard Feynman diagram expansion provides the necessary infor- 
mation on the structure (position and shape) of the peaks in the matrix elements. In 
general, the main source of variation of the amplitude in a given region of the phase 
space can be associated with one (or more) of the propagators in one (or more) dia- 
grams becoming large. Typically, this happens either when collinear-divergent limits 
are approached or when virtual massive particles, with a small width, are close to 
their mass-shell. 

At present, the available codes for the automated generation of phase space are 
based on identifying singularities in Feynman diagrams || [7|. As we will show in 
the following, we use a similar approach for generating the phase-space mappings. 
However, our method differs from the previous ones in that we propose to exploit 
the information gathered from the Feynman diagrams to reorganize and simplify the 
way we perform the integration over the various channels. 

In this paper we present a new multi-channel integration method, which we will 
call Single-Diagram-Enhanced, and its implementation in the new multi-purpose 
event generator MadEvent which is powered by a new version of MadGraph ||. Given 
a user process, MadGraph automatically generates the amplitudes for all the relevant 
subprocesses and produces the mappings for the integration over the phase space. 
This process-dependent information is passed to MadEvent, and a stand-alone code 
is produced that allows the user to calculate cross sections and to obtain unweighted 
events. Once the events have been generated (and all the necessary information 
stored 1 ) they may be passed to any shower Monte Carlo program (such as HERWIG [ID| 



or PYTHIA |0, O) where partons are perturbatevely evolved through the emission 
of QCD radiation, and eventually turned into physical states (hadronization) . 

This paper is organized as follows. In the next section we present our method 
for mapping the phase space and integrating the squared amplitude. In Section [|, 
we provide some examples of the results that can be obtained with MadEvent. The 
conclusions and a brief discussion of the outlook are addressed in the last section. 



addition to the usual kinematic quantities, such as masses and momenta, the spin, the color 
connection, and the flavour of each of the final-state partons may be needed for each event. This 
information is stored in the "Les Houches" format |B|. 
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2. Single-Diagram-Enhanced multi-channel integration 



In this section we present our approach for integrating a generic squared amplitude 
over the phase space of the final-state particles. We first review the usual multi- 



channel integration method, following closely Ref. ||13| , since this serves to establish 
the notation and allows us to point out the issues relevant to the subsequent discus- 
sion. 

— * 

Let us consider /(<&) the function to be integrated over the set of phase-space 
variables $. Given N the number of particles in the final state, $ is of dimension 
d = 3N — 4. In order to account for the complex structure of / several mappings 
of the Monte Carlo variables x G [0, l] d are introduced and each of them generates 
a different probability density, = 1, ...,n, where n is the total number of 

channels. 2 The optimization of the relative importance of the different channels is 
achieved by introducing a set of weights ati, Y^i=i a i = ^i a i > 0> an d defining a 
total probability density g($>) = Yli=i a i9i{®)-> which combines the weights of each 
channel and is still correctly normalized. The weight of a phase space point $ can 
be written as = /($)/p($) and the result of the Monte Carlo integration and 
its variance are: 



1=1 d$g($)w($) = / (2.1) 

V(a) = [ d$g($)w($) 2 = [ (2.2) 
J J g($) 

As made explicit in the notation, the result of the integration does not depend on the 
choice of the (which can even differ iteration by iteration), while the variance does. 
The aim of the multi-channel approach is to find the best set of a« which minimizes 



the variance V. Several strategies can be followed |T3[. Here, we just comment on a 
few important aspects of this approach. 

In the first place we note that for each point in the phase space, independently 
from the actual channel chosen, the determination of g(<&) implies the computation 
of the contributions from all channels. As long as the number of channels is not very 
large and/or the computational burden for calculating g($) is much smaller than that 
for /($), this might not be a serious problem. However, in general, this procedure is 
inefficient. During the integration, the relative weight of each channel changes due to 
the optimization of the a{. In particular, even if some of the channels become less and 
less important compared to others, the time spent to calculate their contributions 
remains the same. Another point of weakness is that, even if the evolution of each 
of the can be chosen to be independent from the others, in practice the a« are 



2 Each probability density <?j(3>) is related to the mapping through $ = G t 1 (x), where G,($) is 
the primitive of gi(<&). 
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correlated. Predicting the behavior and/or the stability of an integration with a large 
number of channels is very difficult and jeopardizes the robustness of the solution. 3 
We argue that there is a natural way to overcome the above limitations and to 
calculate the integral of a squared amplitude over the phase space in an efficient way. 
The idea is very simple and is based on the following observation. Assume that the 
function to be integrated could be written in terms of a basis of n functions /, such 
that: 

n 

/ = X)/< with /;>0, Vi, (2.3) 
1=1 

and such that the peak structure of each can be efficiently mapped by a single 
channel gi. Then, the integration of / reduces to: 

/ = Jd^fm = ^Jd$g l ^)^ = ^I lJ (2.4) 

i.e., to a sum of n independent integrations. For a generic integration problem, the 
identification of such a basis might be too difficult and the above method be simply 
not viable. However, in the case at hand, a natural decomposition arises from the 
physical content of the process. We propose to use the following basis: 

fi — \—\ I a i 2 |-^tot| ) (2-5) 

where is the amplitude corresponding to a single Feynman diagram and A tot = 
J2i A is the total amplitude. It is clear that each fi in Eq. (|2.5|) satisfies the condi- 
tions (|2.3|) and form a "complete" basis. 4 More importantly, the peak structure of 
each ft is the same as of the single squared amplitude \Ai\ 2 . Finding the suitable 
mapping is therefore straightforward, since it can be derived from the propaga- 
tor structure of the corresponding Feynman diagram. We will refer to the above 
algorithm as Single-Diagram-Enhanced multi-channel integration. 

Decomposing the integration of the amplitude into n independent integrations 
has immediate advantages. First, in contrast to the standard multi-channel integra- 
tion, the evaluation of the weight in one channel does not require the computation of 
the others. This entails that, from the statistics point of view, in the Single-Diagram- 
Enhanced integration method the complexity of the computation does not increase 
with the number of channels. 

Second, the Single-Diagram-Enhanced integration is parallel in nature. Integra- 
tions in each channel can be performed by using different resources and the results 

3 This aspect can in practice hamper the possibility of using many channels if the mapping gi 
are obtained by a stratified Monte Carlo method. 



4 Eq. is analogous to the procedure suggested by Odagiri 14 to extract the leading color 

behaviour of QCD amplitudes. 
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combined only at the very end. Considering the latest developments in this field, 
such as the availability of PC farms at low cost, we reckon this as an alluring property 
of this approach. 

A third useful aspect of our method is that, given a target accuracy of the 
integration, it is trivial to reweight the channels so that those whose contributions 
to the total result is small (large) are evaluated with a smaller (greater) number of 
Monte Carlo points. One simply introduces the analogues of the a, as in a standard 
multi-channel method, defining them as on = h/I. This realizes in practice the 
desirable property that the computational time spent on a single channel should be 
proportional to its relevance to the final result. 

Another interesting property of the decomposition ( |2.5|) is its intrinsic modular- 
ity. The number of the diagrams to be included in the decomposition is completely 
arbitrary and dictated only by convenience. In general, the number of channels will 
be less than the number of Feynman diagrams. The first important example is when 
k identical particles are present in the final state and n can be easily reduced by a 
factor k\. Another case is when some diagrams do not display any peaked behav- 
ior and are dropped from Eq. (|2.5| ). Alternately, some diagrams can be efficiently 
grouped together if the form of the resulting is predicted by general arguments, 
such as radiation coherence. 



3. MadEvent at work: some examples 

In this section, we present some results for processes relevant at hadron colliders, 
which have been obtained with MadEvent. The purpose is to show the versatility 
and the potential of our approach and also to provide benchmark cross sections 
for the user. Detailed information on how to use the MadGraph/MadEvent package 
together with a web interface for generating the code for a given process is available 
at [http : //madgraph . physics . uiuc . edu| . In order to facilitate the comparison with 



other available codes we have chosen a very simple set of cuts and conventions, which 
are summarized in Table [p. In particular, the cuts in transverse energy and AR are 
left quite loose so that regions of the phase space where the matrix elements peak 
are not left out and the efficiency of our integration technique can be fully tested. 5 
In Table [2] we present the cross sections involving the production one or more 
vector bosons in association with jets, both at the Tevatron and LHC. Apart from 
the interest they might have on their own, these processes contribute to some of 
the most important backgrounds for top-quark studies (PU+jets, Wbb+jets) , for the 
Higgs searches (W /+ W / ~+jets, Wbb+jets, ZZ+jets) and for SUSY particles searches 
(Z+jets, with Z — > vv). In the single vector boson production, final states are 



5 In fact, at the LHC such cuts are not sufficient to ensure a good behaviour of the perturbative 
series for some of the cross sections shown, as it is apparent by comparing the processes with n jets 
with the corresponding ones with (n + 1) jets. 
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leptons and so that the vector bosons are not contrained to be on-shell. This entails, 
for instance, that cross sections involving a Z boson, correctly include the diagrams 
where 7* — > £ + £~ and their interference with Z* — > £ + £r . 

For simplicity, in the case of double vector boson production, we present the 
results where the vector bosons are produced on shell and are not decayed. Nonethe- 
less, we stress that MadEvent correctly handles also the decays, where non-trivial 
phase-space mappings have to be introduced. An example is provided by W + W~ 
production where contributions from diagrams where ^*/Z — > W + W~ — > £ + £~vv 
appear, which demand particular care. In this case, there are different regions in the 
phase space which can give large contributions, coming from the Z and one of the 
W's, or both Ws being on shell. We have ensured that the procedure used for the 
mapping of the phase space automatically accounts for all the possibilities, so that 
no approximation is introduced. 

When possible, we have compared our results with the ones obtained by running 
available codes (such as VECBOS @, Wbbgen |H| and MCFM Q) and found good 
agreement. We plan to conduct a detailed set of comparisons with the recently 
released Alpgen [|Tj| in the future. 

Table ^| contains the cross sections for the production one or more heavy-quark 
pairs in association with jets, both at the Tevatron and LHC In the case of tibb 
production we also give the cross sections at order o% m a 2 s . These last ones consti- 
tute the irreducible backgrounds in the search of the Higgs when it is produced in 
association with a top-quark pair. Also, the cross section for two bottom-quark pairs 
in association with jets are listed. In the case of two bottom-quark pairs we have 
also performed a detailed comparison with the results of Ref. 0] and found good 
agreement. 

In Table |] cross sections are given for the production of an Higgs boson in 
association with a top-quark pair, with jets ("W-boson fusion") and in association 
with a W boson, both at the Tevatron and LHC. 



4. Conclusions and Outlook 

In this paper we have introduced a new integration method, called Single-Diagram- 
Enhanced multi-channel integration, and its implementation in the multi-purpose 
event generator MadEvent. 

We have discussed how the problem of integrating over the phase space ampli- 
tudes that involve a rather large number of particles in the final state can be reduced 
to the integration of amplitudes associated to single Feynman diagrams. The lat- 
ter are much simpler tasks and can be tackled by many of the techniques already 
present in the literature. Our method best applies to the automated computation of 
cross sections (and, more importantly, to the generation of unweighted events) for 
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colliders 


Tevatron: pp @ 2 TeV 
LHC: pp @ 14 TeV 


QCD and EW couplings and parameters 


a s (M z ) = 0.1185 

a e m(M z ) = 1/128.9, sin 2 ^ = 0.2312 
CKM matrix: diagonal 


p.d.f. 


CTEQ5L (parametric) 


scales 


Hf = Hr = M z 


vector boson masses and widths 


m z = 91.188 GeV, T z = 2.495 GeV 
m w = 79.96 GeV, T w = 2.06 GeV 


heavy-quark (pole) masses 


m t = 174.3 GeV, m b = 4.7 GeV 


Higgs parameters and couplings 6 


M H = 120 GeV, T H = 3.7 • 10" 3 GeV 
y t = 0.97, y b = 0.017 


kinematic cuts 


E T > 10 GeV, \r]\ < 2.5, AR > 0.4, 
for all the final-state particles 



Table 1: Couplings and conventions used in the calculation of the benchmark cross sec- 
tions. 



processes where many different parton subprocess contribute, such as vector boson 
production in association with jets at hadron colliders and many others. 

We have presented some applications of MadEvent, a multi-purpose event gen- 
erator based on the above ideas, to studies relevant at the present and forthcoming 
colliders. All the necessary steps to generate a set of unweighted events are per- 
formed automatically. The user only inputs the process and the desired integrated 
luminosity or the total number of unweighted events. 

Although in principle there is no limitation to the number of particles that can 
be considered in the final state, in practice this still depends on the complexity of 
the process itself and on the total number of Feynman diagrams involved. Currently, 
the package is limited to ten thousands diagrams per subprocess. So, for example, 
W+5 jets is feasible but close to the present limit. 

Several improvements that will extend the ability of MadGraph to deal with a 
larger number of QCD partons are under study. One necessary step is the "factor- 
ization" of the amplitude in order to reduce the complexity of the calculation from a 
factorial growth to an exponential one, as already achieved for non-QCD amplitudes 
in Refs. fH 0, |2| and extended to the QCD amplitudes in Refs. @, |2(| [Ij§. In 



the latter case, a reorganization of the color structure is also needed and the sum 
over the colors of the external QCD partons can be efficiently performed by a Monte 
Carlo integration. Work in this direction is in progress |23| . 



6 In fact Th and yi, do not enter the results presented here, but we prefer to leave them as a 
reference for the user. 
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Table 2: Benchmark cross sections for single and double vector boson production in 
association with jets at hadron colliders. The decays into leptons, W — > ev e and "y*/Z — > 
e + e~, are included in the single vector boson cross sections and cuts, as described in 
the text, are applied to all the final state particles (jets and leptons). In the case of 
WW, WZ, ZZ production, the vector bosons are not decayed and are produced on shell. 
Statistical errors for all processes are at the percent level. 
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Table 3: Benchmark cross sections for heavy-quark pair production in association with 
jets (light quarks and gluons) at hadron colliders. The heavy quarks are produced on shell 
and are not decayed. In tibb and bbbb production, contributions from the exchange of a 
virtual Higgs have been excluded. Cuts, as described in the text, are applied to all the 
final state particles. Statistical errors for all processes are at the percent level. 
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Table 4: Benchmark cross sections for Higgs production in association with a top-quark 
pair, with jets (EW processes) and in association with a W boson, at hadron colliders. The 
heavy quarks and the Higgs and the W are produced on shell and are not decayed. Cuts, 
as described in the text, are applied to all the final state particles. Statistical errors for all 
processes are at the percent level. 
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